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Abstract 

This study concerns a situation when measurements of the nonresonant cross-section of nuclear 
reactions appear highly dependent on the environment in which the particles interact. An appealing 
example discussed in the paper is the interaction of a deuteron beam with a target of deuterated 
metal Ta. In these experiments, the reaction cross section for d(d,p)t was shown to be orders of 
magnitude greater than what the conventional model predicts for the low-energy particles. In this 
paper we take into account the influence of quantum effects due to the Heisenberg uncertainty 
principle for particles in a non-ideal medium elastically interacting with the medium particles. 
In order to calculate the nuclear reaction rate in the non-ideal environment we apply both the 
Monte Carlo technique and approximate analytical calculation of the Feynman diagram using 
nonrelativistic kinetic Green's functions in the medium which correspond to the generalized energy 
and momentum distribution functions of interacting particles. We show a possibility to reduce the 
12-fold integral corresponding to this diagram to a fivefold integral. This can significantly speed 
up the computation and control accuracy. Our calculations show that quantum effects significantly 
influence reaction rates such as p +'^Be, '^He +^He, p +'^Li, and ^^C +^^C. The new reaction rates 
may be much higher than the classical ones for the interior of the Sun and supernova stars. The 
possibility to observe the theoretical predictions under laboratory conditions is discussed. 
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I. INTRODUCTION 



The rates of non-resonant nuclear reactions including fusion processes are determined by 
kinetic energies of the interacting particles in their center-of-mass system as well as by their 
distribution in energy and momentum. For moderate temperatures the main contribution 
to the fusion process is expected from particles with energies several times larger than the 
plasma temperature. 

However, it is well known that in dense environments the quantum uncertainty in the 
energy of particles associated with their frequent collisions leads to disruption of the unam- 



biguous relationship between the energy and momentum of particles [lH3|. This results in 
the appearance of power distributions in the momentum distribution function of particles in 
dense media. It is particularly interesting to study how these effects in a nonideal plasma 
contribute to the rates of fusion reactions at moderate plasma temperatures of a few electron 
volts and densities of about one gram per cubic centimeter 

The influence of quantum effects on the equilibrium momentum distribution was investi- 
gated by Wigner and others , who found the amendment to the Maxwellian distribution, 
proportional to the square of the Planck constant. Their result was related to the noncom- 
mutativity of the operators of kinetic and potential energy of interaction and suggested an 
effective "temperature" increasing with the density of particles. 

The first results on the power dependence of the equilibrium distribution function (as 
for real Bose-particles were obtained by Bogolyubov in his famous work on superfiuid weakly 
nonideal Bose-gas (see also Q). Similarly, one can obtain the momentum distribution for 
electrons in the superconducting phase at momenta larger than the Fermi momentum. Later, 
the result of the power law distribution function for momenta larger than the Fermi value 

n 

was obtained by Belyakov [9|] for electrons interacting with short-range impurities. In 1960 
Vosko et al [lO\ found the distribution function of electrons above the Fermi momentum, 
taking into account the exchange interaction of electrons. That function decreased beyond 
the Fermi jump as p~^. 

In 1966 Galitsky and Yakimets jll| showed that the equilibrium momentum distribution 
of particles acquires a power correction term to the Maxwellian function due to quantum 
effects. This correction is valid at large momenta that exceed the thermal or the Fermi 
momentum. For the Coulomb interaction potential, the inverse power dependence on the 



momentum was equal to eight. In 1975 Kimball proved that in a Coulomb system the 
asymptotic momentum distribution also includes the eighth degree, regardless of the plasma 
temperature 



121 ] ■ In the collision of particles of one species in theory there is an additional 
factor proportional to the correlation function at zero distance between them. This factor 
may be of the order of 1/2 for electrons colliding with each other due to the Pauli exclusion 
principle. For heavy particles such as hydrogen or deuterium nuclei this factor is small in 
the strongly coupled plasma. 
In 
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151 ] the authors proposed a simple model using the Lorentz gas concept in which 
a light particle is scattered by a heavy impurity particle. It was demonstrated numerically 
that there must be a significant deviation of the distribution function from the Maxwellian 
form as the result of the quantum corrections. This fact leads to the power law damping 
of the distribution function, the exponent is also equal to eight for the screened Coulomb 
interaction. It has been suggested that this effect should influence the reaction rate for 
these particles. In particular, it should give the nonexponential temperature dependencies 
of reaction rate constants for inelastic processes. 
In the literature 
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171 ] there was some criticism, connected with the use of the asymp- 
totical presentation of a single-particle distribution function for calculation of reaction rates, 
including fusion rates. The problem is that in reality we must use the product of distribu- 
tion functions over momentum in the laboratory frame for both reacting particles and due 
to power low tails, not Maxwellian ones, it is hard to perform analytical integration of the 
reaction cross-section, depending on particles relative momentum in their center of mass. 

In this study we simulated the thermonuclear fusion reactions taking into account the 
impact of this mechanism on the distribution function under conditions that might be re- 
alized in a contemporary experiment. The rates of some reactions were calculated. It was 
shown that in general the reaction rate is determined by the diagram shown in Fig. 1 for the 
single-particle nonrelativistic kinetic Green's function or generalized distribution of energy 
and momentum, which corresponds to the escape process of species a and is reduced, in 
general, to a ten-fold integral. In other words, in the original expression, one should not 
perform a simple averaging of the reaction cross-section, depending mainly on the relative 
momentum of particles (rather than energy) for quantum single-particle momentum distri- 
bution function. It is also shown that under certain conditions this integral is reduced to 
a five-fold integral and in the model case to a three-fold integral. In some cases the last 




FIG. 1. Reaction rate diagram for a + 6 — ?■ c + d 

integral can be calculated explicitly and contains contributions from the power type momen- 
tum distributions for each of the reacting particles in addition to the classical Maxwellian 
terms. Each stage of reduction and the corresponding simplification was verified numeri- 
cally without the use of any simplifications. Some attempts have been made to calculate 
several fusion reactions appropriate for conditions of the solar plasma in the deep interior. 
A reasonable accuracy of the simplified explicit estimates for the rate constants was shown 
for these conditions. 

From our calculations we predict in important circumstances a significant increase in the 
rate of neutron yield. In particular, an increase in the rate of d + d reactions may in fact be 
observable in laboratory testing. We also note possible changes in the reaction rate constant 
in the plasma of astrophysical objects. 

For some reactions the influence of the power corrections was shown to be very significant 
at relatively low temperatures and high densities of weakly nonideal plasma. 

II. CALCULATION OF REACTION RATES IN NONIDEAL PLASMAS 

The state and properties of a system are found from the generalized distribution function 
F{E,p) where E and p are the energy and momentum of the particles. This function 
should be defined in a factorized form. Further we use kinetic energy to substitute for the 
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momentum so the generalized function can be written as 



F{E,e) = n{E)5^{E-e) (1) 

where n{E) are the occupation numbers and 5'-){E — e) is the spectral function, describing 
the dependence of the generalized distribution function on energy and momentum. Equation 
(II]) is the most general representation of the nonrelativistic kinetic Green's function {2!, 
without any assumptions, in which the width of the Lorentzian spectral function ^7 is the 
imaginary part of retarded mass operator of the particle in a medium and energy shift is 
the real part of retarded mass operator. Consequently, the reaction rate follows from the 
integration which is formally a twelve-fold construction written as (see Fig. [1]) 

S = ^^^2 J ^P"- Jq J ^P'> J J dq 

(E,) {l±n{E, + Qa-u)) 
X 5-fb (Eb - 6b, Eb) n (Eb) il±n{Eb + co + Qb)) (2) 

X Hb {Eb + UJ + Qb- ep,+q, Ep^+q) Iff . 

Here the subscript indices a and b designate the reacting species. Ea and pa are, respectively, 
the energy and the momentum of particles a, Ea describes the kinetic energy, fiab is the 
reduced mass, u is the interaction energy, h is Planck's constant, and Qa is the energy 
released in the fusion reaction. The dependencies of the occupation numbers and "±" are 
subject to proper statistics. Here "— " refers to fermions and "+" to bosons. The actual 
irreducible dimension of this integral is equal to 10. 

For the nonideal plasma, the distribution function contains the Lorentzian which deter- 
mines the spectral dependence on the kinetic energy: 

Sl{E-E,E) = imZ!L_ (3) 

^ ^ {E-E-A{E,E)f + j^ ^ ' 

The scattering linewidth in the Lotentz gas model 7(£',£:) is found from 

^{E,E) = hNaV (4) 

where is the number density of the scatterers, a is the scattering cross section, and V is 
the collision velocity, determined by energy E. 



In Ref. [18|] it was noted that, for the gaseous medium approximation, when only binary 
colhsions can be considered, the reaction amphtude is just a function of the momenta before 
and after their reaction. In our further transformations we use the value of the reaction 
cross section which depends on the energy in the center-of-mass system. The ratio between 
the amplitude and the cross section has the form 

Iff = ^-/(^.) (5) 
where p is the momentum of the reacting species a and b in the center-of-mass system: 

^ mbPa - niaPb \pf 

p = ; = o — (6) 

The reaction cross section is the function of kinetic energy e„ in the center-of-mass system 



and may be written as in [l9|, i.e. 



^/(^p) = exp {-27r?7(£p)} (7) 



where r] = ZiZ2e^ /hv is the Sommerfeld factor. It is conventional to use a different form of 
the cross-section for a non-resonant fusion reaction: 



a,(^,) = ^exp (8) 



Sp y y Sp 



where the Gamow parameter is found from 



Ea = ^^"^"^:f'^''' = 4,.a,^Z!ZiRy (9) 
n me 



with Ry = 771^6^ /2h . With a very good accuracy we can approximate it as Ry = 100yUabZ^Z| 



[keV]. At the same time, the factor S{Sp) is weakly dependent on energy Sp. 

In a dense medium with account for the effects of degeneracy the scattering amplitude for 
nuclear fusion may depend on the total energy j2o|. This leads to corrections proportional 
to plasma concentration multiplied by the cube of the elastic scattering amplitude, which 
are low within the gas approximation. 

The reaction rate found from ([2]) is a very general definition formulated for this model. 
The modeling procedure is reduced to calculation of the tenfold integral which is a very 
complicated task. Under conditions when plasma becomes ideal, for example for smaller 
densities, we may write 7 {E, s) and ^7 {E — s, s) can be reduced to the delta function. 



For reactions with energy release {Qa > 0), in case of nondegenerate plasma the pop- 
ulation numbers are small and n{E) can be neglected as compared to unity. With this 

simplification the reaction rate obtained from the general equation is the following 
1 r 

^ = 1 YT^ / / dpaha {Ea - Sa, £a) n {Ea) 

X ^ dEb J dpbS'jbiEb- £b,£b)n{Eb) (10) 

X j duj j dq {Ea + Qa-UJ - £pa-q: ^Paq) {Eb + UJ + Qb - £pb+q: ^Pi+q) Iff ■ 

In this approximation we also neglect the suppression of transmission for Fermi particles 
or amplification (condensation) for Bosc particles. This is a valid assumption because of 
small populations in both cases. We further assume that the lincwidths of the energy and 
the kinetic energy distribution profiles, determined by 6j'^ and ^7^, are small enough to 
approximate the contours with the Dirac delta functions. With this assumption in equation 
we come to 

^ = 9 \^ / ^^0. I dpaU {Sa) S'Ja {Ea - £a, £a) 

infj,lf,{2nh) Jo J 
X ^ dEbJ dpbU (sb) 5-fb {Eb - Sb, £b) \p\ o {^p) (11) 

X J duj J dqS{Ea + Qa-uj - £p„_g) S{Eb + uj + Qb- £pj+g) 
The inner integration over uj in this equation can be performed as following: 

Ig^ J dq j dujd {Ea + Qa- <^ - £pa-q) 5{Eb + Qb + u: - Sp^+g) 

= J dq6{Ea + Qa + Eb + Qb- ep^-q - Sp^+g) 
Now it is easy to obtain the ratio for the delta function arguments: 

+ 2,n, + 2mb 

1 1 \ ^ ( ,^ ^ Pb 



= + — + — ]-2q cos {paq) - - — cos {pbq) 

Let us reduce it to the perfect square by introducing the energy in the ccntcr-of-mass system: 

P^ {Q—pf_ {Q^pf /..x 
^pa-q + ^Pb+q ^ £a + £b- 7, ^ 7> — ea + Sb - ep-\ (14j 

so that the integral Ig now takes the form 

Ig^ J dqS {Ea + Qa + Eb + Qb- Sp^-g - Sp^+q) \ (15) 



= j dqS (Ea + Eb + Qa + Qb- £pa 



\q-p\ 

{q-pf^ 



'^fj'ab J \q-p\ 



Let us note that variable Eg as well as momentum p is determined by the momenta of colliding 
particles and does not depend on q. Changing to the variable s = q — pwe get 

\2 \ 1 



Iq = J ds 6 (Ea + Eb + Qa + Qb - Spa - ^Vb + " 
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2/iaf,/ \q-p\ 

/ s^ds 5{Ea + Ef, + Qa + Qb-£p^-Sp,+ep-^]— (16) 

"'O V 2j2abJ \S\ 

POO 

= AnHab / des6 {Ea + Eb + Qa + Qb- £pa - ^Pt + ^P - ^s) = '^T^l^ah- 
JO 

The integration in the latter equation was carried out with the use of the variable Eg = 

\s\'^/2flab- 

Now we substitute the result into equation (fTTj) . For the relative velocity of the colliding 
particles we need to make an account for Vab = Va — Vb = |p|//Xa6- As the result we come to 

r°° dEa f , ^ , ^ , ^ , /""^ dEb 



S = I / dPa n {Ea) {Ea - Ea) f 



X 



/\p\ 
dpb n {Eb) 5^ {Eb - Eb) fXab cr {Ep) Air^ab (17) 

r°^dEan{Ea) f.^r^r^ . dEbn{Eb) f,^r.r^ a t/ / a ^ ^ 
= ^^^^Y J ^ ^ " " Jo {2nhf 1 ~ ■ 

If we use the apparent ratio for the number densities Na and Nb this equation can be 
transformed to the standard form, i.e. 

S = {VabCr) NaNb (18) 

Numerical simulation of a fusion reaction in the light of the above models and the calculation 
of the reaction rates with ([2]) and fll7p as well as comparison of the results confirmed the 
correctness of the approximations used in our approach. Table [T] shows the constants for the 
model reaction (aVab) for values of particles concentration and temperature T = 2 eV in 
the region where quantum effects are important. The reaction constants (crVab) obtained with 
the numerical model in which we used equation (2) is designated as (Vabc) -10 and the results 
obtained with equation (ITTl) are denoted as (\4bO")_5. The detailed simulation algorithm in 
the Monte Carlo method is described in the appendix. These results were compared with 
the analytical estimations based on (73) which is obtained later in this paper. Calculations 
were performed for the astrophysical factor S (0) = 241 keV barns and the Gamow energy 
was determined in accordance with the equation (j9]) as Eq = 100 keV. 
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N = 10^1 cm^ 


N = 10^2 cm =^ 


N = 10^=^ cm^ 


Analytic 


1.5 X 10-29 


1.5 X 10-28 


1.5 X 10-2^^ 


{Vabcr) -5 


1.3 X 10-29 


2.0 X 10-28 


3.2 X 10-27 


{Vaba) -10 


0.99 X 10-29 


2.1 X 10-28 


2.7 X 10-27 



TABLE I. 

III. KINETIC ENERGY DISTRIBUTION FUNCTION WITH QUANTUM COR- 
RECTIONS 



For the nondegenerate plasma the occupation numbers should have Maxwellian distribu- 
tion over energy at high temperatures: 

= (-tI;) (19) 



V^ikTf ' \ kT 

At the same time the energy distribution function could be derived by integrating the gen- 
eralized distribution function as shown here: 

f{e) = f{E,e)dE= n (E) 5^{E - e) dE (20) 

= '-rn (E) '^llE4 ^ dE. 

71 Jo ^ ' {E-e-A{E,e)f + j{E,ef 

For the values of the kinetic energy in the range of e < Eq, where the threshold energy is 

significantly greater than the linewidth of the Lorentz function, i.e. 7 << Eq, the Lorentzian 

can be represented by (5-function with a high accuracy. Therefore, for this energy range the 

kinetic energy distribution function is f{e) = n{e). If the kinetic energy is e > Eq, then in 

order to calculate the integral one should split the integration interval in Eq. fj20|) into two 

regions, i.e. the region of small energies of the order of the plasma temperature Eq : 3kT...5kT 

and the rest of the interval. In the region oi E > Eq the Lorentzian can be approximated 

by the delta function. In the region of small E with account for 7 << kT in the asymptotic 

limit of e >> kT the denominator of the Lorentzian is approximately e'^: 

Wo \E-e-A{E,e)f + ^{E,ef 

+ r n (E) , dE (21) 

Jeo ' {E-e-A{E,e)f + j{E,ef ^ ' 
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= U {e) + /o (£) . 

The classic expression for the distribution function is represented by the term 

2 



In the range of small E it is necessary to use the Coulomb cross section: 



^t{e,) = (23) 

and the expression for width of the Lorentz function which is conditioned by scattering of 
particle a on plasma particles I: 

-ial = llNiGt[ep)Val (24) 

In the center-of-mass system we write 

£p = /^az ( — + — -24/— — cos(p;,pz) I (25) 
\ma mi \ niami ) 

If we consider the Lorentz gas approximation, i.e. nia <S m^, then 

E £ 

^a/rUa » Ei/mi, — (26) 

fJ'al rUa 



where 



I - |2 
\Pa\ 

2ma 



(27) 



or Va,^J^ (28) 

V f^ab \ \ ma 



The cross section for Coulomb scattering of particles a by particles I: 



Oai {e) = = ^^-^ (29) 



. J5 (30) 



where = (m^ + m/) f« 1 for « m^. 
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FIG. 2. The distribution function of the kinetic energy in the Lorentz gas for temperature 0.01 
KeV (a), obtained from numerical simulation of the generalized distribution function (b) and the 
Maxwellian distribution (c). 

In a multicomponent medium the line width is determined by the sum of the contributions 
from different sorts of scatterers (strictly speaking, the sum over / must take into account 
the particles different from species a, according to the Lorentz model): 



-,„ = i:-,„, = _^^_^iV,Z,-5- = — (31) 



Here we introduced the following notation: 



2 



^ai = Y.NiZl-i (32) 



In order to calculate the "tail" of the distribution function, we need to substitute the line 
width from Eq. ( l30l) into integral ( 12T]) and write 

- /^xp (-^ j (33) 



^ {kT) 



The influence of the correction to the distribution function is significant only in the asymp- 
totic region of its argument which is the kinetic energy. It should be noted that at low values 
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FIG. 3. The share of the quantum correction to the full distribution function in the approximation 
of the Lorentz gas for temperature 0.01 KeV (a), the same quantity for the generalized distribution 
function, obtained from numerical simulation (b), the share of the Maxwellian part to the full 
distribution function (c). 
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FIG. 4. The same as in Figure [2] for temperature 1.0 KeV. 



of the kinetic energy e < Eq the distribution function is determined by the classical expres- 
sion, i.e., the Maxwellian function ( l22l) . For the temperature of 10 eV and the concentration 
of interacting particles (deuterium ions) = 10^'^ cm~^ the value of the energy threshold 
parameter can be defined as Eq = 5kT. In this case, taking into account equation fl30l) one 
can male the following estimation: 



TiNaV 2TihNe'^ 2kT 1 



{kTf 



12 



m Eq 



5 X 10"^ 



(34) 




log(s) 



FIG. 5. Curves (a)-(c) the same as in Figure [2] (temperature 0.01 KeV), curves (d)-(f) respectively 
the same functions of temperature 1.0 KeV. 

Figure 2 shows the contribution of the quantum-tail correction in the approximation of the 
Lorentz gas in the full distribution function. Figure |2] shows the Maxwellian distribution 
functions with a correction that takes into account the quantum effects in the approximation 
of the Lorentz gas and obtained in the analytical form fl?I]) . as well as the distribution func- 
tion of the kinetic energy, obtained by numerical integration of the generalized distribution 
function. 

In Figure[3]this contribution is presented as the ratio fqt{E) / f{E) in equation (pT]) . where 
- curve (a), the same value for the distribution of kinetic energy, obtained by numerical 
integration of the generalized distribution function - curve (b), as well as the contribution 
of the Maxwell distribution in the full distribution function kinetic energy, the ratio of 
the formula fl2T]) - curve (c). The simulation was performed under conditions of a shock 
compression of deuterium in the titanium matrix with the density of atoms of 5.7 x 10^^ 
cm~^ at a temperature of 0.01 keV. As one can see the energy is less than 0.1 keV the total 
distribution function is determined by the Maxwellian term. At energies higher than 0.2 
keV the contribution of the quantum correction to the total distribution function becomes 
crucial. Quantum correction obtained in the approximation of the Lorentzian gas as well as 
obtained by numerical integration of the generalized distribution functions, represented by 
curves (a) and (b) respectively. 
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Figure S] shows the same distribution function as in Fig. [21 but at a temperature of 
1.0 keV. As one can see from comparison of these figures, with increasing temperature the 
quantum effect comes into play at higher energy. 

Figure O shows the energy distribution for different temperatures: 0.01 and 1.0 keV, 
curves shown are collected from Figs. H] and |H Given that the characteristic energies are of 
different orders, we had to switch to a logarithmic scale. It is clear that in the asymptotic 
region the tails of the distribution function do not depend on temperature and are practically 
identical. From this pictures it may be estimated the percentage of deuterium ions in the 
asymptotic region of distribution function over momentum. This number is much larger 
then correspondent number in Maxwellian distribution, which is considered for calculating 
fusion rate constant in classical approach. 



IV. COMPARISON OF KIMBALL'S APPROACH AND THE LORENTZ GAS 
MODEL 

The equations for the power-tail distribution function of particles momentum in the 
Lorentz model can be compared with the rigorous result obtained with the use of the Kim- 
ball's approach. Considering the repulsive Coulomb interaction between two particles with 
charges Za and Z^, and reduced mass mah from the Schrodinger equation it follows: 



1 o ao ' 

\ Zfiab Tab 

where F does not contain singularities at small r. At short distances one can solve the 
equation (extraction of the singularity at small r) 



where aab is the Bohr radius: 



^ = / 1 + (36) 

V O-ab) 

Using this solution one can come to the distribution function by calculating the Fourier trans- 
form of the product of the solutions with singularities at the particle species a approaching 
particles / and m. In the final expression points / and m are to tend to each other. As a 
result, we obtain the asymptotic momentum distribution function in the following form 

fai ip) = -j^. ^ (38) 
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It is easy to see that these asymptotics in the framework of the Lorenz model and the 
Kimball's approach are in good agreement. This agreement takes place independently of the 
mass of the particles - whether electrons or nuclear particles - and whether or not they are 

degenerate. 

It is interesting to note that, in contrast to the Lorenz model, the Kimball's formulas 
are valid for an arbitrary mass ratio of the colliding particles. For the scattering of particles 
of species a on one another we can get, using Kimball's method, the following: 

/..(P) = ^^«..(0) (39) 

In a nonideal plasma these contributions are small in proportion to the ion-ion correlation 
function. Contribution to the quantum asymptotic of the distribution function for protons 
and other ions at expense of their scattering on neutrals and electrons is small due to small 
cross sections and the square of the reduced mass. 

One can compare the asymptotics of the distribution function obtained by Kimball and 
within the Lorentz model for the interaction potential, which has a singularity near zero 
distance between the particles: 

(n < 3) (40) 
The Schrodinger equation with such a potential is written as 

^ ^'v^, + ^V = ^ (41) 



and has a solution 



After calculating the square of the Fourier transform of this solution we obtain 

^p-;^!^ (43) 

At the same time, from the expression for the quantum correction to the Maxwellian distri- 
bution we get 



(44) 



where 



atr^k^""-^ (45) 
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Calculating the scattering amplitude on the potential of this type in the Born approximation 
we obtain: 



which is in agreement with ( 143|) . 

In addition to the agreement, as noted above, of the asymptotic expression fl38|) to the 
limit obtained by Kimball for the electron momentum distribution function due to their 
interaction with ions, with the results of the Lorentz model, this treatment now allows us to 
generalize the result for the distribution function of heavy particles such as deuterons inter- 
acting with heavy ions, as soon as they satisfy the Schrodinger equation and the Coulomb 
law at distances of the order of the Bohr radius fl37|) . Note that the result obtained by 
Kimball does not depend on the ratio of masses of interacting particles, since the interaction 
is considered in the center-of-mass frame. Thus, we can conclude that the power law for the 
distribution function in the asymptotic region holds not only where the Lorentz model is 
applicable to interacting particles, but also for arbitrary particles. For example, it holds also 
for the interaction of deuteron ions with hydrogen ions. This asymptotic behavior does not 
depend on the temperature, so the Kimball's theory for electron distribution function and 
Lorentz gas model considered here coincides in spite of great differences of physical objects, 
i.e., metals near zero degree in the first case and dense plasmas of 1-1000 eV temperature 
in the second. 

It may also be noted that if electrons are localized within their Bohr radius, with an 
uncertainty in their momentum and kinetic energy of the order of 13.6 eV, then for deuterons, 
similarly localized to a Bohr radius fl37|) . the scale of uncertainty of their kinetic energy will 
be of the order of 50 keV, which is on the scale of energy necessary for fusion to occur. 

V. REACTION RATE FOR THE MODEL DISTRIBUTION FUNCTIONS 

The reaction rate for the case of the nondegenerate plasma at sufficiently high temperature 
is determined by averaging the interaction frequency of the particles over their distribution 
functions. The equation to calculate the reaction frequency is reduced to finding 




(46) 




(47) 
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Here Vab is the relative velocity of particles a and b, Sp is the energy of the particles in the 
center-of-mass system. 

The relation between the energy of the particles in the center of mass coordinate system 
and in the laboratory system has the form: 



gp = /iafe f — + — - 2 J cos {pa.Pb) I (48) 



In the laboratory system we have 

Vab = Va-Vb 



^ = 72. 



— + — - 2.p^ COS {pa.Pb) (49) 



^ab \ TUa rUb ^ maVflb 

Using equation ( l47l) for calculating the reaction frequency might be easier if one reduces the 
integral's dimension by performing the angular integration: 

/■oo p /• 1 

{aV)=2{2'Kf dpaplfiSa) jdpbplfiSb) j_^d cos 9 abVabCTfiSp) (50) 

Using the same approximation for the distribution function ( |2T]) . i.e. marking out the 
power asymptotics, we can rearrange the equation to form convenient for standard numeric 
integration. The reaction frequency is then split into four terms, which correspond to the 
terms of the distribution function in the region of small energies and in the asymptotic 
region. 



(aV) = 2(27r) / dpaPa (/o (Sa) + fqt {Sa)) / dpbPb ifo (eb) + fqt {Sb)) / d COS OabVabCTfiep) 

Jo Jo J-1 

roo roo rl 

= 2(27r)^y^ dpaPafoi^a) dpbPb f oi^b) J ^d COS 6abVab (T f (Ep) 

+ 2(27r) dpaplfo{ea) dpbPbfqt{eb) J ^d cos OabVab a f{ep) (51) 

roo roo rl 

+ 2(27r)M dpaPafqtiSa) dpbPb foi^b) d COS 6 abVab Ci f{ep) 

Jo Jo J-1 

roo roo rl 

+ 2(27r)M dpaPafqtiSa) dpbPb fqt{£b) d COS 6 abVab CT f {Sp) 

i/ <^ iJ — 1 

This equation for the reaction rate constant can be used in calculations of fusion reaction 
rates. Along with equations ([2]) and f|T7|) this expression is the next order approximation 
suitable for numeric modeling of fusion reactions. 

We change the variables of integration from the momenta to energies and substitute 
expression ( |33l) for the asymptotic distribution function. As a result we obtain: 
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+ 



J "^aV^exp|-— W -jj^ J d cos 9 abVab (T f [Sp 
^ 



~ — .3/2 " / ~m f d^by/£bexp(-^\ I d cos 9 abVab cr f (ep) 



H , / / ^772 / d cos 9 abVab(^f{ep) (52) 

.Jm~ffrb Jeo el'^ Jeo el' ^-i 

Further calculations of constants for the reaction with the use of this model are denoted 
as (Vabc) _3. The last formula takes into account that the part of the distribution function 
conditioned by the quantum effects is fqt{£) = in the range of low kinetic energies. 

Equation (l52l) can be analyzed to estimate the causes and conditions under which the 
asymptotic region of the energy distribution functions gives a determining contribution to 
the rate of fusion reactions. The integrand in the first term of f l52|) is a series of factors 
that are notably different from zero in different domains of their arguments. The values of 
the Maxwellian distribution function are ~ 1 for the energies not too much higher than the 
temperature. In the asymptotic region these functions decrease exponentially. The fusion 
cross section is exponentially small at low energies and reaches its maximum for the energies 
of the order of hundreds of keV, as, for example, happens in the synthesis of the deuteron. 
Therefore, at low temperatures the first term in (!52|) becomes small. Note that it is the one 
which determines the reaction rate in the model currently accepted for computation. 

The second and third terms are equal in calculating the rates of reactions of identical 
particles, such as deuterons. Despite the small factor preceding the integrals these terms 
may exceed the first term for sufficiently large concentrations of the particles. Let us note 
that according to equation ( l32l) T^ab ~ N, therefore the contribution of the asymptotic 
distribution function increases with increasing density. In these terms the main contribution 
to the reaction rate is given by the different definition regions of the integrands. For example, 
for the second term the integration of the Maxwellian function over energy Sa is substantial 
for the values of the argument only slightly exceeding the temperature. In order for the value 
function - the reaction cross section was not negligible the domain of its argument should 
be taken in the order of hundreds of keV. Hence, the integral over Sb gives a significant 
contribution to the asymptotic region, where the function decays as a power law. Thus, for 
estimations it is possible to assume that the reduced energy in the center of mass Ep does 
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not depend on Sa- The integral to this argument is evaluated as: 

, / rf^av^exp ~ 1 (53) 

In order to estimate the last or the fourth term in (152|) we can calculate the ratio of its value 
to the value of the second term. Given the recent relation obtained we get: 

2he''Zl'Lai r dea _ 2he''Z^,J:ai 2/5 



el^^ ^5/2- 



(54) 



This ratio is obtained for the values T^at ~ ~ 10 cm and Eq = 5kT = 5 ■ lOeV. Note 
that the latter estimate is consistent with the relation (13^ 



VI. REACTION RATE IN A LORENTZ GAS 

Let us perform calculation of reaction rates conditioned by different energy ranges (the 
argument of the distribution function). The reaction frequency for particles of species a and 
b in the Lorentzian gas approximation, i.e. motionless target particles, is reduced to the 
following calculation: 

roc 

{ctV) = f{e)Vaba{e,)^ede (55) 
Substituting the distribution function, we obtain: 

{^V) = / (/o {e) + U {e)) Va {e) ^ede (56) 

Calculating the terms separately, we get: 



T 



'id iw) 

The reaction rate determined by the asymptotic part of the distribution function is found 
from 

W^)qt = J fqt i^a) Vab 0"/ {Sp) ^adSa = 

y/rna Jo ]l IXab I V J ^a^P 

iiiS(0)rexp|-..m^<fe„. 

^af^ab Jo [ V J y/^^a 



Ahe'ZiT. 
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When calculating the rate of fusion of identical particles, such as dd reaction, in the center- 
of-mass system for particles of equal masses and equal energies we have: 



^ab = /^a6 ( — + — - 2 J — — cos {Pa.Pb) ) = '^l^ab— (1 " COS {Pa,Pb)) 

\ma nib 



= 2- 



rria + rua rUa 
Finally, we come to 



1^ IJlamb J rUa 

COS {Pa:Pb)) = £a (1 - COS {Pa.Pb)) ~ ^a- 



/rUa 



5(0), 



ma/ 2 Jo 



exp < —7r\ 



de 



(60) 



(61) 



4v^;ie4z2E„i5(0) 



m„ 



exp < — TT^ 



de 



AV2r,c'ZjE,,SiO) 2-5! 8 • 5! v^r/c^ZfE^^ 5(0) 



The intermediate integral used in this equation has been calculated as following: 

de 



n^Elm„ 



= 2 



exp 



-7r\ 



IEg 



Jo 



X e 



-dx 



x-^ 



Jo 



dx 



(62) 



= 2^ 



9q;^ V q; 



2^fl 



/-I 



/ a=TT^/EQ 

2-5! 



2-5! 



In this calculation we performed the change of variables: 



X — — — ^ dx — 
Ve 



1 7 7 

-as, de = —2— 



(63) 



2^3/2 

In the center-of-mass for particles of identical masses but different energies, such as >> £(, 
we have 



^ab l^ab 

Finally, for such a case we get 



+ --2i 
ma mb ' 



' ^b / -> \ \ 
COS [Pa, Pb) ^ iJ^ab 

ma mb I ma 



(64) 



/m, 



ls(0) / 

a V ^^ab Jo 



f 

Jo 



exp 



IEg \ ^de /fj^y/^ 
£^v^ \maJ 



o4,,2 -72 



exp < —Ti\ 



IEg 



de 



Ahe^/ii^Z-^EiSiO) 2 • 5! _ 8 . 5! 3(0) he^ //^.Z^E. 



(65) 



Given the equal masses of particles we get 

8 • 5! S{0) fie^ mlZl^i _ 2 • 5! S{0) fie^ Zl^i 



7r6 Eg El Ami 



7r6 Eg El ma 



(66) 
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VII. ANALYTICAL ESTIMATES FOR REACTION RATES FOR MODEL DIS- 
TRIBUTIONS WITH QUANTUM EFFECTS 



For the integrals in the second and third terms in (1521) . different ranges of Sa and Sh are 
essential. For the second term the integrand is close to if the energy is greater than the 
temperature, i.e. outside the interval < < 3A;T For variable Eh this integral area is 
significantly wider. In the third term variables Sa and Sb interchange. 

Let us now consider the second term. We assume Sa « £b taking into account the 
significant range of energies of different particles. Then the energy in the center of mass can 
be represented by 



ep = ^iab\— + — - 2J — — cos {pa.Pb) 1 ~ ^^ab— (67) 

\ma ruh \ niamb J nib 

If we substitute this expression into the second term in (152 p it turns out that the rate and 
the cross section do not depend on the angle between the velocity vector and the integral is 
reduced to 

^"^^"^ v^ikTf' v^b Jo ^^"^^"prfcT/io 



— -jj^Vab{ep)(r{ep) (68) 



/ab (£p CT 



^JrfTb KrribJ Jo ej"^ 

2^/2he^ZlT.bi //ioty/' dep_ /^^ex 

y/m^ Knib) Jo V /"ab £p 
Ahe^ZlT^bi fUahV de„ f I E, 



IEg 



••b \mbJ Jo 



^7^— exp<! -tt/ 



TUb XTTib/ Jo ep'\/e; [ V J 

7r6 Eg E'^KmbJ rrib ' 
We can practically perform similar transformations as earlier but here we have 



Z!Ebi = Z^mlY^'-^ (69) 

I ^^bl 



For the third tern in fl52l) we get 



8-5!g(0) he' f fiabY Z'^Jl^i 
7r6 Eg ElKmJ m„ 



(70) 
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where 

AT y2 

Zl^ai = ZlmlY.^ (71) 
Thus, the reaction rate corrected with the quantum tail contribution is equal to 



In the sums of the last formula the summation index / should stand for all particles of 
the medium except for species a ox h respectively. It follows from ( |39l) and the following 
comments. Because of the small correlation function of identical particles, the quantum 
corrections corresponding to such scattering are also small. 

Let us note that the reaction rate constant determined for the tails of the distribution 
function and represented as (1721) does not depend explicitly on temperature. As noted above, 
this is due to the fact that in the asymptotic region the tails of the distribution function 
do not depend on temperature. At the same time it should be noted that temperature 
defines the ionic composition of the plasma. Therefore, the temperature dependence of the 
distribution function is extremely important and is determined by the terms in brackets in 
( 1721) . As it is shown below, the definition of the ionic composition of plasma is essential for 
temperature range of 5-10 eV. It must be noted, that correct definition of the distribution for 
ionic composition is very complicated problem in nonideal plasmas, so the measurement of 
neutron yield gives an instrument to have information on equation of state in such plasmas. 

Thus, we have given proof of the validity of the approximate approach, which consists 
in averaging the cross sections found using the quantum corrections to the momentum 
distribution functions of the reacting particles. The accuracy of this approach is discussed 
further below where the results of calculations are shown for the full formula ([2]), simplified 
equation ( 1T7|1 . and approximate analytical estimations ( 157|) and ( ITOl) . This eliminates the 
issues raised in the paper 16|] on the validity of such a method. 

From f l72|) we can estimate the rate constant for fusion of deuterons taking 49.6 keV barn 
for an astrophysical factor 5(0) and to 99.9 keV for the Gammov energy Eq. Calculating the 
scattering cross section at Eq as well as estimating the velocity of the particles appearing in 
the reaction rate, we find that the coUisional width of the gamma at the density of scattering 
particles of the order of 3 x 10^'^ cm~^ appears to be about 10-6 eV and the reaction rate 
constant turns out to be of the order of 10-28 cm'^ c^^. 
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le reaction rate is specified to 



13, 
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423l|. The corresponding 



In a coupled plasma the Maxwellian contribution to t 
make an account for screening of the Coulomb potential 
generalized expression is obtained, for example, by replacing the exponential Sommerfeld 
factor in equation ([7]) with the semiclassical tunneling probability through the screened 
potential barrier: 



exp {-27rr]} ^ exp { -h!^ £ ^ _ Hir) - E,dr ) (73) 

where r„ and r/ are classic stopping points and ZaZi^e^ /r — H (r) is the interaction potential 
with allowance for screening effects. It may be noted that the progress in the theory of the 
subbarrier tunneling, connected with the problem of heavy ion fusion, was analyzed in 24 1. 

Great interest was aroused by the experimental studies of reactions involving light nuclides 
such as isotopes of hydrogen. He, and Li. They have shown the exponential enhancement of 



the nuclear reactions cross-sections at low energies |25N30| . The beam of the light nuclides 
ions was focused on the metallic target. The Coulomb screening of electrons had been 
considered in this case. The electrons respond by accumulating around the positive charge 
and therefore partially screen out its positive Coulomb potential. Although the experiments 
have proved the significance of electron screening, a theoretical explanation is still far from 
satisfactory. The high screening potential value arises from the environment of the light 
nuclides in the metallic matrix, but a quantitative explanation is missing. The screening 
effect was far beyond the expected value. 

It must be noted that for the quantum correction the screening effects can be ignored 
because the screening energy is of the order of H{0), and is therefore small compared to the 
Gamow energy Eg, which determines the "tail" contribution to the value of reaction rates. 



VIII. MODELING OF FUSION REACTIONS 



Nuclear fusion reactions occur at an appreciable rate in conditions of hot plasma. This 
is primarily the plasma in stars; in particular, it is the plasma in the solar interior. From 
the models considered in this paper, the quantum effects can also predict the increase in the 
rate of reactions in plasmas of moderate temperatures but high densities. By monitoring 
the synthesis reaction, this increase might be observed for suitable parameters. In order 
to determine these parameters, we calculated the fusion rates for different conditions. It 
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N = 3.311025cm-=^ 
1.3362 KeV 


N = 9.5754311024cm^ 
0.6892 KeV 


N = 7.4051022cm^ 
0.1828 KeV 


Classic rate 


2.9810-^5 




0.7910-51 


Quantum correction 


2 7710-49 


0.8 • 10-^9 


0.62 • 10-51 


Full rate (analitic) 


2.98 • 10""^^ 


6.13 • 10^"^"^ 


1.41 • 10-51 


{^Vab) -3 


3.27 • 10-^5 


7.45 • lO-"^"^ 


1.28 • 10-51 


{(yVah) -5 


3.08 • lO-''^ 


6.78 • lO"'''^ 


1.43 • 10-51 


{aVab) -10 


3.35 • lO-''^ 


6.48 • 10-"^^ 


1.61 • 10-51 



TABLE II. 



is also interesting to do some revision of processes in the interiors of stars, to re-evaluate 
the contribution of various processes to the release of fusion energy, and ultimately to the 
evolution of stars. 

Table II shows the constants, which were calculated for the fusion reaction p + p — >■ 
D + + at three points along the solar trajectory, from the central part to the periphery. 
In the table we show the constants calculated numerically: (Vafe)_3, (V^b)_5, (Kf,)_10 as well as 
analytical estimations of constant for the reaction, which were calculated using equation (58) 
(Classic rate), the quantum correction of constant for the reaction, which were calculated 
with equation (73) (Quantum correction), and the sum of these two values (Full rate). As 
one can see the calculations performed for different models are in a satisfactory agreement 
with each other. 

Another conclusion to be drawn from the results presented in this table is that the 
influence of the quantum effects is notable for this reaction only in a distant region from 
the Sun center. But in this region the fusion rate is much lower than in the central region, 
so the influence of the quantum corrections brings no effect on the energy balance for this 
reaction and synthesis of deuterons in the solar interior. Our simulations were carried out 



using the data on the astrophysical factor from |31 |. 

An example of the reaction with the rate greatly influenced by the quantum corrections 
is the reaction of the hydrogen cycle: ^He +^ He — ?■ 2p +^ He. Table III shows the rate 
constants calculated for the same points of the solar trajectory. In calculations of these rates 



the value of the astrophysical factor was taken from 
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N = 3.31 lO^^cm^ 
1.3362 KeV 


N = 9.575431 • lO^^cm 
0.6892 KeV 


N = 7.405 • lO^^cm =^ 
0.1828 KeV 


Classic rate 


6.84 • 10-34 


5.65 • 10-39 


1.86 • 10-53 


Quantum correction 


2.34 • 10-27 


5.37 • 10-28 


4.14-10-30 


Full rate (analitic) 


2.34 • 10-27 


5.37 • 10-28 


4.14-10-30 


{oVab) -3 


2.38 • 10-27 


5.42 • 10-28 


4.15 - 10-30 


{aVab) -5 


1.34 • lO-^'^ 


3.03 • 10-28 


2.24 - 10-30 



TABLE III. 



Similarly, it could be shown that the rates of many fusion reactions, which occur in the 
solar interior, such asp +^ Be, ^He +^ He, etc., as well as of reactions such as C + C in the 
depths of supernovae, become much larger than their classic values if the quantum effects 
are taken into account. These results require further consideration and beyond the scope of 
this work. 

However, these predictions, in principle, might be verified in laboratory experiments using 
a dense plasma with moderate temperatures around one electron volt and above. This can 
be achieved in explosive experiments like the ones, for example, that study the equation of 
state of strongly coupled plasma |33l-l35|. If we take deuterium, compressed to a pressure of 
the order of megabars, we can make the following prediction: in pure deuterium, in which 
there is a noticeable degree of dissociation and ionization due to pressure ionization, the 
reaction rate will be very small in these conditions because of the factor g{0) in ( 139|) . which 
in turn is small in a strongly coupled plasma. If deuterium is diluted with an extraneous gas, 
such as a different isotope of hydrogen, helium, etc., then there will be terms in the D + D 
reaction rate due to scattering on a buffer gas, which can lead to an observable neutron yield 
of about 10® to 10^° neutrons per pulse of about 1 microsecond length. The result strongly 
depends on the ion composition and possibility of the plasma phase transition in strongly 
coupled plasmas. 

As another example, we present calculations for a mixture of deuterium and xenon, which 
can create a shock wave with a speed exceeding 5 km/ c and create a plasma with temperature 
of 5 — 10 eV in the reflected wave. Calculations of the plasma parameters for these conditions 
were kindly provided by V.K. Gryaznov. 
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T, eV 


9.27 


9.88 


1.06 


1.22 


N(D+) 


5.45 • lO^i 


5.54 • 10^1 


5.63 • 102i 


5.74 • 102i 


N(Xe+) 


6.35 • 1020 


2.79 • 10^0 


1.46 • 1020 


6.24 • 10^9 


N(Xe+2) 


5.53 • 10^1 


5.19 • 10^1 


4.61 • 1021 


3.31 • 1021 


N(Xe+=') 


2.26 • 10^1 


2.68 • 10^1 


3.08 • 102i 


3.76 • 102i 


N(Xe+4) 


8.82 • 10^9 


1.40 • 10^0 


2.19 • 1020 


4.79 • 1020 


{aV) 


3.5 • 10-29 


3.68 • 10-29 


3.85 • 10-29 


4.13 • 10-29 


(aF)N(D+)2 


1.04 • 10^5 


1.13 • 10^5 


1.22 • 10^5 


1.36 • lOis 



TABLE IV. 



Table HV] shows the results of the calculation of constant fusion reaction: D+D 3He+n 
for different spatial points of the plasma in conditions of shock compression of a mixture of 
deuterium and xenon. The gas mixture was D : Xe = 50 : 50 in volume ratio. The initial 
pressure was 25 bars. The table shows the equilibria concentration of various components 
of plasma for these conditions: N(D"'"), N(Xe"'") etc, reaction constants {(jV) and the rates 
of the reaction (ar)N(D+)2. The latter ones have been obtained analytically by taking into 
account the quantum corrections. 

As this table shows, the temperature of the plasma in such conditions is more than 2 
orders of magnitude lower than the temperature of the Sun plasma. The plasma density is 
sufficiently high. Under the conditions of shock compression of Xenon is 2 to 3-times ionized 
as a result of pressure ionization 32| . The calculations used data on the astrophysical factor 
of the reaction of 36 1. 

For the conditions of the shock experiments it would be important to estimate the re- 



laxation time of the distribution function. In [3|, it was suggested, that the only binary 
elastic scattering must be taken into account to study relaxation of the nonequilibrium dis- 
tribution function. It is the Coulomb type collision for the plasma in a shock wave. The 
cross section for these collisions decreases with increasing energy. So the elastic scattering 
frequency for the kinetic energy is of the order of the Gamow energy. When the fusion 



reaction is realized distinctly, it may be estimated as: 



2vre^Z2z2 

Eh ' 



2K 



G 



(74) 



Then the upper level of the relaxation time is r ~ [N {cftV)) ^ with the data from Table HVl 
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In this case it may be estimated as about 1 ns to 10 ns. In [37|] the period of relaxation of the 
nonequilibrium distribution function was calculated more precisely. The model, used in this 
paper, permits an analysis of the relaxation process for all velocities, including formation 
of the distribution function tail. In this paper for the various dependencies of the elastic 
scattering cross section on energy it was obtained that the time of the equilibrium setting 
is larger, up to factor 6, than that of the elastic collisions. This estimation was obtained for 
the energy equal to the kinetic temperature. The calculated value of relaxation is less than 
1 ns. 

Thus, for estimating the fusion reaction, it can be seen that the collisions in the plasma 
of the shock experiments are effective in bringing the distribution function to equilibrium 
during the hydrodynamical process with the characteristic time of < 1 fis. 

The fusion rate constants for the above conditions calculated with no account taken for 
the quantum effects are by about 20 orders of magnitude smaller than the values listed in 
the table. The reaction rate for different points of the plasma vary by ~30%, as it is seen 
from the last line of the table. The lifetime of such a plasma is of ~1 /xs, which gives the 
neutron yield of ~1.2 x 10^ cm~^. 



IX. CONCLUSIONS 

In this paper we analyzed the influence of quantum effects on the rate of fusion reactions. 
As a result of frequent collisions of particles in a dense plasma there disappears the com- 
plete correspondence between the total and kinetic energy of the particle, the generalized 
distribution function, thus, depends on both the total and the kinetic energy of the particle. 
The momentum distribution function has the power dependence on the kinetic energy in the 
asymptotic region. Carrying out the averaging over the distribution function to calculate, 
for example, the fusion reaction rate leads to a notable increase in reaction rate as compared 
to the calculations using the Maxwellian distribution function. We created numerical models 
for various conditions of the reactions and carried out calculations in a wide range of plasma 
parameters. The use of the approximate analytical estimates obtained under the averaging 
procedure has been validated. 

In addition, the ranges of parameters where the most pronounced quantum effects are ex- 
pected are shown. Two experiments are identified that might demonstrate quantum effects: 
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one, by comparing DD reaction rates in densely compressed cold deuterium plasma with or 
without a buffer gas; and, two, by arranging for a shock wave in mixtures of Deuterium and 
Xenon. 

The presence of quantum tails is also evident in the vibrational kinetics of low-temperature 
plasma. It increases the rate of V-T relaxation, as shown in [3|. This result of the theory 
is in a very good agreement with the experimental data under normal conditions. The 
theory also predicts the reduction of the induction time for ignition of hydrogen-oxygen and 
hydrogen-air mixtures at pressures above 5 atm and at temperatures below 1000 K. 
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Appendix A: Monte Carlo integrations 

Calculation of the approximate value of the integral using the Monte Carlo method is one 
of the few methods for calculating the quadrature in this problem, given the high multiplicity 
of integral. Optimization of the calculations in the framework of this method is the choice 
of determining the probability distribution of sites of integration of the quadrature formula. 
On one hand it minimizes the dispersion of the mean value while on the other hand it 
minimizes the possibility of a fairly simple and efficient simulation of random vectors, which 
determine the quadrature grid of the numerical integration. 

Minimum dispersion in the calculation of the integral is achieved when selecting a random 



distribution function of several variables, proportional to the integrand [38|]. Given this con- 
dition, the probability density of the random vector in a multidimensional space was defined 
as the product of functions, which were probability densities for individual components of 
the random vector. In the Monte-Carlo method the integration variables are the very ran- 
dom variables. Functions that are selected as probability densities constitute a significant 
part of the integrand function provided in the multiplicative form. For example, we consider 
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in detail the computation of the integral, which reduces the calculation of the reaction rate 
for the non-degenerate distribution function ( !T7|) : 



rOO nOO roc pOC pi 

Sab= / dEa dsay/e^ I dEb deb^/el dcos{pa,Pb) (Al) 

JO JO JO JO J-1 

X n{Ea)S-fa {Ea - Ea, Ea) n{Eb)S-fb {Eb - Eb, Eb) Vab CT (Ep) . 

In accordance with the foregoing notes, in integration over the variables Ea, Eb we use 
n{Ea), n{Eb) as probability density functions, and over the variables Ea, Eb - functions of 
the spectral particle characteristic S'jaiEa — Ea,Ea) and SjbiEb — Eb,Eb)- 

In accordance with the technique used in this paper we will decompose the factors into 
terms corresponding to the probability density and a factor, the average value of which is to 
be calculated. The following expressions describe for the probabilities of random variables, 



i.e. the variables of integration flAip with the selected probability density: 

^3a {Ea) = a;a f^'^ dE^u (E,) , < E„ < oo (A2) 

JO 

S2a{ea)=Ol2a dE aS'Ja (Ea - E a, Ea) , < E a < OO (A3) 

JO 

Ssb (Eb) = asb r' dEbU (Eb) , < Eb < oo (A4) 

JO 

S2b (Eb) = «2b / dEbS'yb {Eb - Eb, Eb) ,0 < Eb < oo (A5) 

JO 

(x) = aia J dx, — 1 < X < 1. (A6) 

Here (z) is the probability of a random variable to be in the range with the upper 
limit of z, the factors are the normalizing factors that ensure the implementation of the 
normalization condition for the probability of the variable denoted by the index ip. 

Distribution of random variables with a given probability density is performed using the 
standard methods. It is necessary to find a solution for the system of nonlinear equations: 

(z) = u^, (A7) 

where index ip runs over all the values corresponding to different variables of integration, 
are random variables uniformly distributed over the range [0,1]. Variables Ea, Eb, x can 
be determined from one of the equations in system (]A2|) . (]A4|) . and (1A6|) . The selected 
probability densities for variables Ea, Eb contain other integration variables in addition to 
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their "own" variables. Therefore, during simulation of these random variables, we have to 
solve the system of two nonlinear equations (IA3|) and (IA5I) applying the method of iterations. 
The value of integral (lAip is found as the average value of the expression: 

1 ^ 

•S'afe = — XI V^V^yab(j{£p), (A8) 

with variables Ea and Eb as well as values of Vat and (y{Ep) contained in this sum determined 
by the above-described drawing of random variables. 
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